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Abstract 

The  use  of  gradient  based  optimization  algorithms  in  inverse  design  is  well  established  as 
a  practical  approach  to  aerodynamic  design.  A  typical  procedure  uses  a  simulation  scheme  to 
evaluate  the  objective  function  (from  the  approximate  states)  and  its  gradient,  then  passes  this 
information  to  an  optimization  algorithm.  Once  the  simulation  scheme  (CFD  flow  solver)  has 
been  selected  and  used  to  provide  approximate  function  evaluations,  there  are  several  possible 
approaches  to  the  problem  of  computing  gradients.  One  popular  method  is  to  differentiate 
the  simulation  scheme  and  compute  design  sensitivities  that  are  then  used  to  obtain  gradients. 
Although  this  black-box  approach  has  many  advantages  in  shape  optimization  problems,  one 
must  compute  mesh  sensitivities  in  order  to  compute  the  design  sensitivity. 

In  this  paper,  we  present  an  alternative  approach  using  the  PDE  sensitivity  equation  to 
develop  algorithms  for  computing  gradients.  This  approach  has  the  advantage  that  mesh  sensi¬ 
tivities  need  not  be  computed.  Moreover,  when  it  is  possible  to  use  the  CFD  scheme  for  both 
the  forward  problem  and  the  sensitivity  equation,  then  there  are  computational  advantages.  An 
apparent  disadvantage  of  this  approach  is  that  it  does  not  always  produce  consistent  deriva¬ 
tives.  However,  for  a  proper  combination  of  discretization  schemes,  one  can  show  asymptotic 
consistency  under  mesh  refinement,  which  is  often  suflicient  to  guarantee  convergence  of  the 
optimal  design  algorithm.  In  particular,  we  show  that  when  asymptotically  consistent  schemes 
are  combined  with  a  trust-region  optimization  algorithm,  the  resulting  optimal  design  method 
converges.  We  denote  this  approach  as  the  sensitivity  equation  method. 

The  sensitivity  equation  method  is  presented,  convergence  results  are  given  and  the  approach 
is  illustrated  on  two  optimal  design  problems  involving  shocks. 
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0078  and  F49620'93- 1-0280,  and  by  the  National  Aeronautics  and  Space  Administration  imder  NASA  Contract  No. 
NASl-19480  while  the  second  author  was  a  visiting  scientist  at  the  Institute  for  Computer  Applications  in  Science 
and  Engineering  (ICASE),  NASA  Langley  Research  Center,  Hampton,  VA  23681-0001. 
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1  Introduction 


Optimal  design  problems  consist  of  selecting  design  parameters  for  a  system  in  order 
to  optimize  a  given  design  objective,  usually  constrained  to  satisfy  a  partial  differential 
equation.  In  many  of  these  problems,  design  parameters  describe  the  shape  of  an 
object.  Examples  of  these  shape  optimization  problems  include  drag  reduction  [21], 
[22],  weight  minimization  [14],  optimal  sensor /actuator  placement  [6],  airfoil  design 
[16],  [17],  [18],  [19]  and  the  design  of  wind  tunnel  elements  [15]. 

Traditionally,  approximate  solutions  of  these  problems  are  found  by  “cut  and  try” 
methods,  combining  a  designer’s  engineering  experience  with  repeated  experimental 
testing.  This  is  often  expensive,  motivating  computational  methods  which  compute 
the  optimal  design  directly.  These  methods  require  defining  an  objective  function 
and  an  appropriate  PDE  model  of  the  states  of  the  system.  A  comparison  of  several 
optimal  design  methods  may  be  found  in  [13]. 

Many  popular  approaches  couple  a  gradient-based  optimization  algorithm  with 
function  evaluations  provided  by  a  proven  simulation  scheme.  One  of  the  disadvan¬ 
tages  of  these  approaches  is  the  expense  of  computing  the  gradient.  Using  finite 
differences  is  often  too  costly,  even  if  appropriate  step  sizes  can  be  found  and  the  sim¬ 
ulation  scheme  can  take  advantage  of  “nearby”  solutions  (as  is  the  case  with  iterative 
solvers  for  nonlinear  equations). 

Two  strategies  for  alleviating  the  computational  expense  of  gradient  evaluations 
are  adjoint  variables  [17]  and  design  sensitivities  [14].  Adjoint  methods  are  advanta¬ 
geous  when  either  the  problem  is  self-adjoint  or  there  are  a  large  number  of  design 
parameters.  However,  when  there  are  relatively  few  design  parameters,  using  design 
sensitivities,  quantities  which  describe  the  influence  of  the  design  parameters  on  the 
states  of  the  system,  is  an  attractive  alternative.  In  addition  to  efficient  gradient 
computations,  they  can  be  used  in  some  problems  to  construct  an  effective  update  of 
the  approximate  Hessian  for  quasi-Newton  optimization  algorithms,  e.g.  [10]. 

A  standard  approach  often  used  to  compute  the  design  sensitivities  is  based  on 
(implicitly)  differentiating  the  simulation  scheme  (for  the  states)  with  respect  to  the 
design  variables.  Using  the  chain  rule  to  carry  out  this  calculation,  results  in  an 
efficient  numerical  scheme  for  the  sensitivities.  The  efficiency  arises  from  reusing 
many  of  the  quantities  computed  in  the  simulation  scheme.  In  fact,  the  “inversion” 
of  the  system  matrix  (i.e.  the  matrix  factorization)  can  often  be  reused. 
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A  disadvantage  of  this  approach  is  that  for  shape  optimization  problems,  the  dis¬ 
cretization  is  parameter  dependent.  Thus,  derivatives  of  the  discretization  (mesh 
sensitivities)  are  required  for  each  shape  parameter.  Depending  on  the  simulation 
scheme  used  for  the  states,  determining  the  discretization  can  require  the  solution  of 
a  partial  differential  equation  (as  is  the  case  for  finite  difference  solutions  of  viscous 
flow  problems  [26]).  This  requires  a  strategy  for  computing  the  mesh  sensitivities 
[20],  or  for  computing  an  approximation  to  them  [24],  [25]. 

Another  approach  to  finding  design  sensitivities  relies  on  approximating  the  par¬ 
tial  differential  equation,  known  as  the  sensitivity  equation.  This  equation  is  obtained 
by  implicitly  differentiating  the  (infinite  dimensional)  state  equation  with  respect  to 
each  design  parameter.  As  shown  in  [2],  using  the  same  numerical  scheme  to  ap¬ 
proximate  the  sensitivity  equation  which  is  used  to  approximate  the  states,  leads  to 
an  efficient  scheme  with  similar  computational  advantages  as  the  design  sensitivity 
approach  described  above.  Furthermore,  since  the  discretization  is  applied  directly 
to  the  sensitivity  equations,  no  sensitivity  of  the  mesh  is  required.  The  sensitivity 
equation  is  always  linear  in  the  design  sensitivity,  even  if  the  state  equation  is  non¬ 
linear.  Since  there  is  no  requirement  to  use  the  same  numerical  scheme,  it  is  possible 
to  gain  additional  computational  savings  by  using  a  scheme  which  takes  advantage 
of  the  linearity  in  the  sensitivity  equations. 

An  apparent  disadvantage  of  this  approach  is  that  it  does  not  compute  consistent 
derivatives.  In  other  words,  the  sensitivity  equation  approach  does  not  capture  the 
sensitivity  of  the  truncation  errors  in  the  scheme.  Thus,  there  is  a  concern  that 
providing  an  optimization  algorithm  with  an  approximation  of  the  gradient  of  the 
infinite  dimensional  objective  function  instead  of  the  gradient  of  the  approximate 
objective  function  would  cause  the  algorithm  to  fail.  One  might  expect,  however, 
that  if  the  gradients  are  “close  enough”  to  the  true  gradients,  then  the  optimization 
algorithm  should  still  converge.  We  show  that  this  convergence  can  be  established  if 
one  combines  compatible  simulation  and  optimization  schemes. 

Trust-region  optimization  algorithms  are  constructed  to  be  globally  convergent  by 
minimizing  a  model  of  the  objective  function  in  a  region  where  the  model  is  “trusted”. 
This  leads  to  robust  algorithms  capable  of  handling  inaccuracies  in  the  model.  In  fact, 
convergence  results  have  been  given  for  these  algorithms  when  the  model  is  based  on 
inaccurate  gradient  information  [7],  [8].  The  results  hold  provided  the  gradients  satisfy 
a  given  error  condition.  Therefore,  it  is  natural  to  consider  an  optimal  design  method 
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which  couples  a  trust-region  optimization  algorithm  with  gradients  computed  using 
the  sensitivity  equation  approach.  We  denote  this  combined  sensitivity /trust-region 
algorithm  by  the  sensitivity  equation  method  (SEM). 

In  this  work,  we  present  and  analyze  the  sensitivity  equation  method.  The  method 
can  be  applied  to  a  wide  class  of  optimal  design  problems,  including  those  mentioned 
above,  however,  we  focus  on  the  particular  example  of  shape  optimization  of  Euler 
flows  in  order  to  illustrate  the  method.  In  the  next  section,  we  describe  two  design 
problems.  In  Section  3,  we  present  the  sensitivity  equation  method  including  the 
trust-region  algorithm  and  the  use  of  the  sensitivity  equation  to  find  the  design  sensi¬ 
tivities.  Furthermore,  we  compare  various  numerical  approximations  of  the  sensitivity 
equation  with  approaches  based  on  the  discretized  equations.  Section  4  discusses  a 
number  of  convergence  issues  and  includes  a  convergence  theorem  for  the  sensitivity 
equation  method.  In  Section  5,  we  use  a  one  dimensional  duct  design  problem  to  de¬ 
scribe  the  implementation  of  the  sensitivity  equation  approach.  Finally,  we  describe 
the  implementation  and  perform  shape  optimization  for  a  two  dimensional  forebody 
simulator  design  problem  where  the  steady-state  Euler  equations  are  used  to  model 
the  state  variables. 

2  Illustrative  Examples 

We  present  two  optimal  design  problems  below  which  are  used  to  illustrate  the  sen¬ 
sitivity  equation  method.  These  problems  consist  of  determining  shape  parameters 
which  produce  a  solution  to  the  Euler  equations  that  matches  a  desired  flow  “as 
closely  as  possible.”  The  first  problem  is  motivated  by  the  design  of  a  wind  tunnel 
element  in  order  to  produce  a  desired  flow  in  the  test  section.  We  study  a  two  di¬ 
mensional  analogue  of  this  problem.  The  second  problem  consists  of  prescribing  the 
cross-sectional  area  of  a  one  dimensional  duct  to  produce  a  duct  flow  which  matches 
a  desired  flow  profile.  This  problem  was  used  by  Frank  and  Shubin  [13]  in  their  study 
of  optimal  design. 

2.1  Forebody  Simulator  Design  Problem 

The  Arnold  Engineering  Development  Center  (AEDC)  operates  a  free-jet  test  facility 
which  is  used  for  full-scale  testing  of  commercial  and  military  aircraft  engines.  Engines 
are  evaluated  for  performance  and  safety  under  various  free  flight  conditions.  While 
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Figure  1:  Forebody  Simulator  Design  Problem 


this  facility  is  large  enough  to  house  engines,  it  is  not  large  enough  to  house  an  entire 
aircraft  forebody.  Thus,  the  effect  of  the  aircraft  forebody  on  the  engine  inlet  flow 
profile  must  be  simulated.  One  way  of  doing  this  is  to  replace  the  actual  forebody  by 
a  smaller  object,  called  a  forebody  simulator  (FBS).  The  use  of  the  FBS  is  illustrated 
in  Figure  1.  The  FBS  design  problem  is  to  specify  the  shape  of  this  FBS  so  that 
it  produces  an  engine  inlet  flow  profile  which  is  as  close  to  some  desired  profile  as 
possible  [15].  The  desired  profile  can  be  determined  by  performing  either  a  wind 
tunnel  simulation  or  a  computational  simulation  of  a  model  configuration  resembling 
a  test  condition  of  the  aircraft  engine. 

In  order  to  demonstrate  the  applicability  of  the  SEM,  we  consider  a  two  dimen¬ 
sional  analogue  of  this  problem.  This  problem,  depicted  in  Figure  2,  is  to  find  the 
shape  of  the  curve  F,  which  produces  an  outflow  that  matches  the  outflow  generated 
by  the  original  (longer)  forebody  as  closely  as  possible.  The  flow,  Q  (consisting  of  the 
density  p,  the  momentum  pui  -t-  pvj  and  the  sum  of  the  internal  and  kinetic  energy 
E)  is  modeled  using  the  steady  state  Euler  equations. 


Inflow 


Outflow 


Ix)ng  Forebody  (Data  to  be  Matched) 


Shortened  Forebody  (Shape  to  be  Determined) 
Figure  2:  2D  Forebody  Simulator  Design  Problem 


(1) 

(2) 

(3) 


where  7  is  the  ratio  of  specific  heats  (7  =  1.4  for  air).  Given  a  forebody  simulator 
shape  r,  the  flow  (^(r)  is  determined  by  solving  the  Euler  equations  (1)  in  the  test 
cell  domain  fl(r)  subject  to  the  boundary  conditions  (for  supersonic  flow): 


Q 

(u,  u)  •  n 


Qia  at  the  test  cell  inflow, 

0  and 

0  at  the  walls  (no  flow  penetration), 
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(4) 

(5) 

(6) 


where  n  and  i  are  the  normal  and  tangential  vectors  at  the  boundary,  respectively. 
The  set  of  admissible  forebody  simulator  shapes  is 

A={Te  C\a,  6)|  r(a)  =  Fa,  T{b)  =  Tb  and  T(x)  >  \/x  e  (a,  6)}  .  (7) 

A  statement  of  the  design  problem  is  given  below. 

Problem  2.1  (Forebody  Simulator  Design)  Let  Q  be  a  desired  flow  at  the  out¬ 
flow  (engine  inlet), 

S  =  {{x,y)\  X  =  b,  Ti,  <y  <c} .  (8) 

Define  the  objective  function 

J{T)=jjQ(r)-mdS,  (9) 

where  Q(r)  represents  the  solution  of  (1)  with  boundary  conditions  (4)-(6)  in  the  test 
cell  fl(r).  The  forebody  simulator  design  problem  is  to  find  F*  G  >1  such  that 

JiT.)<J{T)  for  all  T  &  A.  (10) 

Closed  form  solutions  to  (1)  with  (4)-(6)  are  available  only  for  special  domains. 
Therefore,  we  consider  approximate  solutions  of  (1)  and  hence  the  approximation  of 
Problem  2.1. 

The  discretization  is  performed  by  selecting  mesh  points  in  the  flow  domain  fl(F) 
where  the  flow  variables  will  be  approximated.  It  is  desirable  to  select  this  mesh 
in  such  a  way  that  the  points  are  more  dense  in  regions  where  flow  gradients  are 
expected  to  be  “large”  (in  order  to  have  more  accurate  differencing)  and  more  coarse 
in  regions  where  the  flow  is  nearly  constant  (in  order  to  save  computer  time).  Other 
issues,  such  as  selecting  points  with  no  sharp  changes  in  density  and  with  sufidcient 
resolution  to  treat  the  boundary  conditions,  make  the  mesh  generation  a  science  in 
and  of  itself  (see  e.g.  [26]). 

Another  constraint  on  the  discretization,  to  simplify  the  implementation  of  a  finite 
difference  scheme,  is  to  use  a  regular  mesh,  i.e.  a  mesh  where  there  exists  a  bijective 
map  taking  the  mesh  points  to  a  lattice  of  points  in  the  computational  space.  For 
example,  suppose  that  A4  is  a  mapping, 

M  :  {x,y)  ^  (11) 

then  derivatives  in  the  physical  space  are  easily  approximated  on  the  lattice  using 
the  chain  rule.  Denoting  the  Jacobian  of  the  mapping  by  Jm-,  the  transformed  Euler 
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equations  become, 


dF  dG  „ 

di  d-q 


(12) 


where 


F  =  C/Q  +  PI 


-1 

M 


’  0  ‘ 

'  0  ‘ 

dr) 

dx 

,  G=VQ  +  PJj^ 

dx 

dr) 

dy 

dy 

_  u  _ 

V 

Q  =  Jm^Q, 


TT  A 

U  =  -^u  +  -;:—v  and 

ox  oy 


dq  dq 

V  =  -^u  +  -^v. 
ax  ay 


(13) 


(14) 


A  standard  finite  difference  scheme,  developed  by  Beam  and  Warming  [1]  is  used 
to  approximate  the  transformed  equations.  The  scheme  introduces  a  time  variable,  t 
as  a  means  of  iterating  an  initial  guess  for  the  solution,  to  a  solution  of  the  steady 
state  equations.  Second  and  fourth  order  artificial  dissipation  terms  are  added  for 
stability,  represented  by  and  respectively.  This  scheme  is  implemented 

in  the  PARC2D  code  [9].  Several  implementation  issues  are  discussed  briefly  below 
which  are  referred  to  in  later  sections.  Readers  interested  in  more  code  details  or  the 
actual  expressions  used  for  and  should  consult  [9]. 

The  difference  scheme  produces  a  system  of  equations  for  the  update  of  the  flow 
variables,  AQ”-  Thus,  the  solution  at  the  nth  iteration,  Q”  is  determined  from 


Qn  ^  Qn-l  ^ 


Sn— 1 


(16) 


The  system  matrix  produced  by  the  approximation  above  is  quite  large  due  to  differ¬ 
encing  in  each  direction.  However,  this  problem  is  circumvented  using  an  approximate 
factorization  into  a  product  of  two  matrices,  each  corresponding  to  differencing  in  one 
of  the  lattice  directions.  The  final  system  has  the  form: 

[/  -F  AtS^A^  -  ^  )  A^Jm]  X 

[/  +  -  V,  Ar,JM]  AO”  = 

-  AtS(F^  -  AtSrfCF 

+  AtV^  (^f )  -  A^  (JmQ^) 

-f  AtVr,  )  -  ^WA,V,)  a,  {JmQ^)  ,  (16) 

where 

-BP  -  -  BG 

Md  B”  =  ^(0“).  (17) 
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The  subscripted  terms  6,  V  and  A  represent  the  central,  backward  and  forward 
difference  operators,  respectively,  in  the  lattice  direction  indicated  by  the  subscript. 
The  converged  solution  is  denoted  by  Q^{x,y)  =  (M(x,y)). 

We  introduce  Bezier  curves  to  parameterize  the  forebody  simulator.  Bezier  poly¬ 
nomials  possess  several  nice  properties  when  used  in  approximations.  The  most  im¬ 
portant  for  the  examples  presented  here  are  the  convex  hull  and  endpoint  interpolation 
properties  (see  e.g.  Farin  [12]).  For  this  problem,  we  consider  a  set  of  two  parameter, 
?  =  (9^5  5^)5  Bezier  curves 

B  =  {r€C>io,i]|  rM  =  (r.w,r,(s;«)),  r,(5;5)>r.,  s€[o,ii,,€]R'’}  (i8) 

where 

r,,(s)  =  aBo, 3(5) -I- O.6B1, 3(5) -h  O.8B2, 3(5)  +  653,3(5),  (19) 

Bj/(5;9)  =  TaBo, 3(5) -1- 3(5) -f- 5^52,3(5) -f  rb53, 3(5),  (20) 

and 

jx’(l-a:)’-’.  (21) 

We  also  assume  a  =  0.5  and  h  =  1.0.  We  can  now  introduce  the  approximate  forebody 
simulator  design  problem. 

Problem  2.2  (Approximate  Forebody  Simulator  Design)  Let  he  de¬ 

sired  flow  measurements  at  S.  JVe  assume  that  the  data  measurements  are  given 
at  the  quadrature  points,  otherwise  interpolation  must  be  used.  Define  the  objective 
function 

Jf(r)  =  Eci  r)-4||],  (22) 

1=1 

where  Q^{xi;  F)  represents  the  approximate  solution  to  (1)  in  the  domain  ^(r)  at  the 
quadrature  point  x,-.  The  approximate  forebody  simulator  design  problem  is  to  find 
r*  €  5  such  that 

J^{T.)<J^{V)  for  all  TeB.  (23) 

Let 

s  =  {(9‘.«“)  s  €  B},  (24) 

then  the  problem  can  be  equivalently  stated  as  finding  {ql,ql)  €  Q  such  that 

foraU  (?',?")  €2.  (25) 
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2.2  Duct  Design  Problem 


This  problem  consists  of  designing  the  cross-sectional  area  of  a  one- dimensional  duct 
such  that,  under  specified  inlet  and  outlet  conditions,  it  produces  a  flow  which  is  as 
close  to  a  desired  transonic  flow  as  possible.  The  governing  conservation  laws  (steady 
state  continuity,  momentum  and  energy  equations)  can  be  reduced  to  a  single  two- 
point  boundary  value  problem  (BVP)  for  the  velocity.  It  was  shown  in  [13]  that  the 
velocity  u,  is  the  solution  of 

Q 

^/(w) =  0,  (26) 

u(0)  =  Uin  and  u(l)  =  Uout, 


where  uin  and  Uout  are  the  velocities  at  the  inlet  and  outlet  of  the  duct,  A  is  the 
cross-sectional  area  of  the  duct, 

f(u)=u+fL^ 


where  H  and  7  are  flow  constants  taken  to  be  1.14  and  1.4,  respectively.  The  Rankine- 
Hugoniot  condition  yields  the  speed  of  sound  as  Ug  =  VH.  Unique  solutions  of  this 
BVP  are  guaranteed  for  monotone  area  functions,  therefore,  cross-sectional  areas.  A, 
are  restricted  to 


A=^Ae  C'(0, 1)  /1(0)  =  A(l)  =  and  ^A[x)  >  0,  Vl  £  (0, 1)  | 


(28) 


for  fixed  inlet  and  outlet  areas  of  Ain  and  Aout*  We  now  describe  the  optimal  design 
problem. 


Problem  2.3  (Duct  Design)  Let  «(•)  €  L^(0,1)  be  a  desired  transonic  flow  proflle 
for  the  duct  and  define  the  objective  function  by 

J(A)=  f  {u{x\  A)  —  u{x)f  dx  (29) 

J  0 

where  u(-;  A)  is  the  solution  to  (26)  corresponding  to  A.  The  optimal  design  problem 
is  to  find  an  A*  G  A  such  that 

JiA,)<J{A)  for  all  AeA.  (30) 

While  the  BVP  has  a  closed  form  solution  [13],  we  consider  approximations  of 
(26)  and  consequently  of  Problem  2.3  in  order  to  study  the  more  general  case.  We 
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begin  by  discretizing  the  duct  length  into  N  cells  (of  length  h  =  j^)  with  centers, 
Xj  —  {j  —  \)h,j  =  1, . . .  ,iV  and  define  to  be  the  average  velocity  in  the  jth  cell, 
i.e. 

=  (31) 


A  system  of  nonlinear  equations  for  (A)  —  found  by  integrating 

(26)  over  each  cell, 

/  fu  (xj  +  a))  -  f  (u  (xj  -  a))  .  , 

^  ^  ^  + « «(^). M^i))  =0.  j  =  l,...,N, 

(32) 

where  it  was  assumed  that  \-§^A  was  nearly  constant  over  each  cell.  An  approxima¬ 
tion  to  is  found  by  replacing  the  fluxes  at  the  cell  edges,  /  (it  {xj  +  |)),  using  the 
cell  center  values  fj  =  /(uf )  and  /j+i  =  /(uJ5.i),  Two  standard  first  order  “Godunov 
type”  methods  are  the  Enquist-Osher  scheme 

/i+i 
/i 

fM 

fj  +  /j+l  “  / (Us) 


/(«(xi  +  | 


jpEO  j 
^j+112  =  < 


jy  ,.N 

3 

N  ^.N 
3 

N 


^3  ’  ^i+1  — 


<  Us  < 


(33) 


<u.<  uf. 


and  the  artificial  viscosity  scheme 


/  +  2  ^  -^i+Va  “  2  +  /i  “  “ (“i+i  “  “i))  5 


(34) 


where  a  has  been  selected  as  1  for  this  study.  These  approximations  were  used  in 
[13],  but  are  included  above  for  completeness. 

We  turn  now  to  the  approximation  of  the  cross-sectional  area  A.  The  space  A 
is  replaced  by  a  subset  of  Bezier  quadratic  polynomials.  The  properties  of  Bezier 
polynomials  allow  us  to  easily  impose  both  the  monotonicity  requirement  and  the 
matching  of  inflow  and  outflow  cross-sectional  areas.  Consider 

B  =  {a  €  (7^(0, 1)  |A(x)  =  AinBo,2(a:)  +  qBi^2{x)  +  Aout^2,2(a:); 

xG(0,l),  ?€  [Ai^,Aout]},  (35) 


where  Bi^r  is  defined  in  (21).  Thus,  B  is  a  one  parameter  set  of  curves  in  A.  We 
restrict  our  optimization  problem  to  this  set  B. 

Our  final  step  in  the  approximation  of  Problem  2.3  is  replacing  the  integral  by  a 
quadrature  rule,  with  the  set  of  quadrature  weights  and  points  {(q,  We  now 

state  the  approximate  design  problem. 
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Problem  2*4  (Approximate  Duct  Design)  Let  represent  data  for  a  de- 

sired  transonic  flow  profile  in  the  duct.  We  assume  that  the  data  and  the  approximate 
solution  are  given  at  the  quadrature  points,  otherwise  interpolation  must  he  used.  De¬ 
fine  the  objective  function 


tzzl 


(36) 


where  u^{A)  is  an  approximate  solution  to  (26)  with  the  cross-sectional  area  A.  The 
approximate  design  problem  is  to  find  an  A^  ^  B  such  that 


J^{A:)<J^{A)  for  all  AeB.  (37) 

Note  that  we  can  identify  any  AeB  with  the  parameter  q  e  Q  =  [Ajn,  Aout]  which 
uniquely  represents  it.  Thus  we  can  equivalently  state  the  problem  as  to  find  q*  e  Q 
such  that 

for  all  q  e  Q.  (38) 


3  Sensitivity  Equation  Method 

3.1  Trust-Region  Algorithms 

We  shall  use  a  trust-region  algorithm  for  the  optimization  loop.  The  reason  for 
selecting  this  type  of  scheme  will  be  clear  when  we  discuss  the  convergence  properties 
in  Section  4.  This  is  a  well  known  algorithm.  However,  we  give  a  brief  description 
below  in  order  to  prepare  for  the  formulation  of  the  sensitivity  equation  method. 

The  quasi-Newton  optimization  algorithm  produces  a  sequence  of  iterates  which 
are  obtained  by  minimizing  a  local  quadratic  model  of  the  objective  function.  This 
model  is  constructed  using  the  evaluation  of  the  objective  function  J^{qk),  its  gra¬ 
dient  (qk)  and  a  secant  approximation  to  its  Hessian,  Hk  at  the  current  iterate 
qk.  The  minimization  of  this  model  produces  the  next  iterate  9fc+i,  i.e. 

mk{qk-^i)  =  min mk{qk  +  Sk)  =  min  {j^{qk)  +  VJ^{qkfsk  +  ^slHkS^  .  (39) 

Thus  the  next  step  is 

tM  =  qk- 

It  is  well  known  that  for  sufficiently  close  initial  guesses  (and  assumptions  on  the 
objective  function),  the  iterates  converge  superlinearly  to  the  minimum,  q^. 


11 


However,  the  initial  guess  may  not  be  in  this  superlinear  region.  Thus  globaliza¬ 
tion  strategies  are  employed  to  bring  the  iterates  into  the  superlinear  region.  It  is 
desirable  to  choose  strategies  which  reduce  to  the  quasi-Newton  algorithm  close  to 
the  minimum.  One  such  strategy  is  a  trust-region  algorithm.  In  this  algorithm,  a 
quantity  8,  known  as  the  trust-region  radius,  is  used  to  measure  the  region  in  which 
the  local  quadratic  model,  m^,  is  “trusted”  as  an  approximation  of  the  actual  ob¬ 
jective  function,  .  Thus,  the  next  iterate,  qk+i,  is  now  found  by  minimizirig  the 
model  in  this  region,  i.e. 


^kiqk+i)  =  „imn  mk(qk  +  Sk).  (40) 

l|Sfcll<*fc 

where  Sk  is  the  trust-region  radius  at  the  kth  iteration. 

A  heuristic  for  changing  the  trust-region  radius  needs  to  be  developed  which  in¬ 
creases  Sk  when  the  model  prediction  is  good  and  decreases  Sk  when  the  model  pre¬ 
diction  is  poor.  One  such  strategy  uses  the  ratio, 

«  =  (41) 

mk{qk)  -  mkiqk+i) 

which  is  the  ratio  of  the  computed  reduction  to  the  reduction  predicted  by  the  model. 
If  this  ratio  is  small  (or  negative),  then  the  model  did  a  poor  job  of  predicting 
and  the  trust-region  is  decreased.  Whereas,  if  the  ratio  is  near  1,  then  the  model  did 
very  well  at  predicting  and  the  trust-region  radius  is  increased. 

We  present  the  resulting  trust-region  algorithm  below. 

Algorithm  3.1  (Trust-Region) 

Select  an  initial  guess  qo  G  Q,  an  initial  trust-region  radius  Sq  and  constants  0  <  rji  < 
7?2  <  1  and  0  <  7i  <  1  <  72.  Compute  J^{qo),  VJ^(go)  and  select  or  initialize  Hq. 
Do  A:  =  0, 1, . . .,  until  “convergence” 

1.  Determine  the  approximate  solution  Sk  to  equation  (40).  We  chose  the  optimally 
constrained  hook-step  method  [11]  to  do  this. 

2.  If  pk  <  Vi,  then  set  Sk+i  €  (0,7i4)  and  qk+i  =  qk,  Jf{qk+i)  =  Jg^{qk), 
VJf  (?jt+i)  =  VJ^^iqk)  and  Hk+i  =  Hk. 

3.  If  <  Pfc  <  7/2,  then  set  G  (0,5fc]  and  ^^4.1  =  qk-{-  Sk-  Compute  Jj^{qk+i), 

and  the  update  Hk+i- 
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4.  If  ri2  <  pk,  then  set  6k+i  G  [4,724]  and  qk+i  =  qk  +  Sk.  Compute  J^{qk+i), 
'^J^iqk+i)  and  the  update  Hk+i. 

Continue 


3.2  Design  Sensitivities 


In  order  to  apply  a  gradient  based  optimization  algorithm,  such  as  the  trust-region 
algorithm  described  above,  we  need  to  consider  methods  for  computing  the  gradient 
of  In  this  discussion,  we  consider  finding  the  gradient  of  (or  a  suitable 

approximation)  with  respect  to  the  single  design  parameter  q.  This  discussion  can  be 
easily  extended  to  find  the  gradient  of  with  respect  to  multiple  design  parameters. 

A  straight  forward  approach  is  to  use  a  finite  difference  approximation,  e.g. 


Unfortunately,  this  approach  may  not  be  practical  for  problems  where  the  approx¬ 
imation  of  the  PDE  is  computationally  expensive,  and  is  overly  complex  in  shape 
optimization  problems  due  to  the  necessity  of  computing  mesh  sensitivities.  One  way 
of  alleviating  the  computational  burden  is  to  use  design  sensitivities,  quantities  which 
describe  the  influence  of  the  design  variables  on  the  flow  variables.  For  example,  we 
can  directly  compute  the  gradient  by  differentiating  (36)  as 

^  JfC?)  =  2  E  Ci  [uf  (?)  -  ti.]  (43) 


The  quantity  is  the  design  sensitivity  for  the  discretized  flow  u^. 

There  are  several  ways  to  compute  this  sensitivity.  As  above,  one  might  use  finite 
differences,  yielding  the  approximation 


dq 


u 


N 


(®»;  ?) 


u 


N 


(xj;  q  H-  Aq)  -  u^jxj;  q) 
Aq 


(44) 


When  the  discretization  is  parameter  dependent,  it  is  easier  to  compute  this  approx¬ 
imation  using. 


d  ^  vF{xi  +  f^M{xi)Aq-,q  +  Aq)-u^{xi\q) 
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in  order  to  avoid  interpolating  back  to  the  unperturbed  mesh.  This  approach  has  the 
advantage  that  it  may  be  possible  to  select  a  step  size  A5  using  error  estimates  for 

u^.  However,  it  is  as  computationally  expensive  as  computing  finite  differences  on 

'tN 
'-'9  ■ 

A  more  efficient  approach  can  be  obtained  by  differentiating  the  simulation  scheme 
used  to  approximate  the  flow  (the  discrete  sensitivity  approach).  For  example,  in 
the  FBS  design  problem,  the  simulation  scheme  (16)  could  be  differentiated  with 
respect  to  9,  leading  to  a  numerical  scheme  for  terms  like  Since  the  chain  rule 

must  be  used  to  carry  this  out,  the  resulting  scheme  for  the  sensitivities  contains 
terms  similar  to  those  found  in  the  simulation  scheme.  Thus,  the  sensitivities  can 
be  computed  efficiently  along  with  the  flow.  A  disadvantage  of  this  approach  is  that 
when  the  discretization  is  parameter  dependent,  as  in  shape  optimization  problems, 
then  derivatives  of  the  discretization  (terms  like  ^Ad)  need  to  be  considered,  see  e.g. 
[20]. 

An  alternative  approach  is  based  on  differentiating  the  original  flow  equation  with 
respect  to  the  design  parameter  and  then  approximating  the  resulting  sensitivity  equa¬ 
tion.  The  result  is  ,  where  the  superscript  N  refers  to  the  approximation 

of  the  flow  equation  and  the  superscript  M  refers  to  the  approximation  of  the  sen¬ 
sitivity  equation.  Since  this  approach  interchanges  the  order  of  differentiation  and 
approximation,  no  mesh  sensitivities  are  required.  Furthermore,  it  has  been  shown 
[2]  that  applying  the  same  approximation  scheme  to  the  sensitivity  equation  leads  to 
similar  computational  advantages  as  the  discrete  approach  described  above.  More¬ 
over,  additional  computational  savings  could  be  obtained  by  applying  a  scheme  which 
takes  advantage  of  the  linearity  of  the  sensitivity  equation.  A  potential  disadvantage 
of  this  approach,  however,  is  that  in  general  ^ 

approximation  scheme  is  used  for  both  the  flow  and  sensitivity  equations. 

However,  if  we  consider  the  gradient  of  the  infinite  dimensional  objective  function, 

^J(9)  =  2^  [u{x-,A)-u{x)]-^u{x;A),  (46) 

then  using  the  sensitivity  equation  approach  provides  an  approximation  of  this  gra¬ 
dient,  i.e. 


Thus,  we  have  reason  to  expect  that  this  approach  could  produce  feasible  gradients 
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for  the  optimization  scheme.  These  two  sensitivity  approaches  are  described  in  detail 
in  later  sections  using  concrete  examples. 

3.3  Sensitivity  Equation  Method 

The  sensitivity  equation  method  couples  a  trust-region  optimization  algorithm  with 
gradient  evaluations  provided  by  approximating  the  sensitivity  equation.  Thus  we 
consider  applying  Algorithm  3.1  with  the  following  quadratic  model, 

M9k+i )  =  II  tpk{qk  +  Sk)  =  II  mm ^  (j^{qk)  +  HkSk^  ,  (48) 

Note  that  we  replace  the  quadratic  model  m*  by  to  emphasize  the  fact  that 
is  approximated  by  gk,  computed  as  {qk)- 

The  intent  is  to  use  the  robustness  of  the  trust-region  optimization  algorithm  to 
compensate  for  the  non-consistent  gradients.  The  result  is  an  optimal  design  method 
which  is  often  more  efficient  and  considerably  easier  to  implement  than  current  meth¬ 
ods.  In  the  sections  below,  we  discuss  convergence  issues  and  describe  the  implemen¬ 
tation  of  this  method. 

4  Convergence  Issues 

Definition  4.1  A  numerical  scheme  is  said  to  produce  consistent  derivatives  with 
respect  to  approximations  N  (for  the  states)  and  M  (for  the  sensitivities)  if 


This  is  exactly  the  case  for  the  discrete  sensitivity  approach,  since  one  actually  defines 
(computes)  to  be 

Definition  4.2  A  numerical  scheme  is  said  to  produce  asymptotically  consistent 
derivatives  with  respect  to  approximations  N  (for  the  states)  and  M  (for  the  sensi¬ 
tivities)  if 

d  /  d  \ 

-  [y/)  ^ 

is  satisfied  as  the  approximations  N  and  M  are  refined. 

We  now  consider  the  convergence  of  the  sensitivity  equation  method.  To  begin 
with,  we  assume  that  the  following  hypotheses  hold. 
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(HI)  For  a  given  qo  in  the  design  space  Q,  let  Qo  he  an  open  convex  subset  containing 
the  level  set  of  at  qo,  i.e. 

Co  =  {qeQ\  J^{q)  <  J^{qo) }  C  Qo  C  Q.  (51) 

(H2)  is  bounded  below 

(H3)  is  Frechet  differentiable  on  Qo 

(H4)  The  Frechet  derivative  of  denoted  by  VJ^,  is  Lipschitz  continuous  on  Qo 
with  Lipschitz  constant  L,  i.e. 

I VJfCs')  -  Vjf (5^)1  <  i  -  «^||  V,',  e  a„.  (52) 

(H5)  The  approximate  gradient,  gu  is  asymptotically  consistent  to  'VJ^{qk). 

(H6)  There  exists  a  constant  ci  G  (0, 1]  such  that 

ci|k||||s,||  <  {-gk,Sk)  <  |b;t|||kfc||  Vfc  =  l,2,...  (53) 

(H7)  There  exist  constants  C2,C3  €  (0,00)  such  that 

-C2  {d,  d)  <  {Hkd,  d)  <cz{d,d)  Vfe  =  1, 2, . . .  (54) 

The  following  discussion  parallels  the  proof  given  in  [7]  which  treats  the  use  of 
trust-region  algorithms  with  inexact  gradient  and  function  values.  This  discussion 
makes  use  of  the  fact  that  we  seek  the  minimum  of  and  have  asymptotically 
consistent  derivatives. 

Lemma  4.1  Under  assumptions  (H6)  and  (H7),  Algorithm  3.1  produces  iterates 
which  satisfy 

i’kiqk)  -  i’k{qk+i)  >  ^Ci||if;fc||  min  1^,  .  (55) 
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Proof  Note  that  since  tpkiqk)  =  J^i'ik), 


'>kk{qk)  -  ^k{qk+\)  =  -  {9k,  Sk)  -  -  {HkSk,Sk) 


(56) 


Now,  let  Sk  =  ll•Sfc||||f^  =  o,*dk,  then  a*  solves 


a  {gk,  dk)  +  {Hkdk,  4)  • 


(57) 


We  can  break  this  up  into  two  cases,  when  {Hkdk,dk)  >  0  and  when  {Hkdk,  dk)  <  0. 
Case  1:  Assume  {Hkdk,dk)  >  0,  then  either 

{9k,  dk) 

*  {Hkdk,dky 


in  which  case 


Mik)  Mlw) 

^  !_(%*>!.>  IjM! 

2  {Hkdk,dk)  ~  2  ^  C3 
using  hypotheses  (H6)  and  (H7),  or 


a^  Sk 


in  which  case 


implies 


Sk  < 


{9k,  dk) 
{Hkdk,  dk) 


M^k)  -  Mqk+i)  =  -4  {9k,  dk)  -  {Hkdk,  dk) 

>  -4  {9k,  dk)  +  ^4  {9k,  dk)  >  ^Ci4|bfc|| 

by  hypothesis  (H6). 

Case  2:  Assume  {Hkdk,  dk)  <  0,  then  a,  =  4-  Therefore 

V’fc(9fc)  -  i’k{qk+\)  =  -4  {9k,  dk)  -  ^Sl  {Hkdk,  dk) 

>  -4  {9k,dk)  >  Cl4||5'jt||  >  ^Ci4|bi:||. 


A 
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Lemma  4.2  Assume  (H7)  holds,  then 


k—i^oo 


imply 


lim 


{sk,gk) 


lim  4  =  0 

k-*oo 

(58) 

=  1. 

(59) 

lkitlllb.il 

Proof  It  was  shown  [11]  that,  if  ||sfc||  =  4,  then  the  solution  to  (48)  is  given  by 
s{fj,k),  where 

^(/^)  =  “  i^k  +  gl)  ^  9k 


and  pk  is  the  unique  real  number  that  satisfies  |k(/i.)||  =  4-  Therefore,  if  4  — >  0, 
then  pk  —>■  oo  (since  Hk  is  bounded,  by  (H7)).  Thus  Sk  —9k^9k-  A 


Lemma  4.3  Let  satisfy  (H3),  (H4)  and  (H7),  then  the  iterates  satisfy 

IMm)  -  *(?wi)]-  -  J"(sHi)]  <  5  (C2  +  i)  IM/,r-(5*  -  . 

(60) 


Proof  Using  the  Cauchy-Schwartz  inequality  and  (H3),  we  obtain 

=  (V  Si,)  +  £  +  A^i)  -  3,}  d\ 

<  (vjf  (»),«)  +  f  ||vjf  (51  +  \S3)  -  VJ«(5()||  llstPA. 

By  the  Lipschitz  hypothesis  (H4), 

(®)  <  {VJ/'(5i),«i)+ /'i||Ast||||si||rfA 

JO 

=  (Vjf(5i).sj)  +  ii|k4|p. 

Thus,  using  (H7), 

[^.(9.)  -  V’fc(9.+i)]  -  [J^{<lk)  -  j/(?.+i)] 

<  -  ^k,3i,)  -  i  {Hus,,  3,)  +  {vj^(qu),3u)  + 

<  -{gu-  Vy/'(56).  Su)  +  i  (C3  +  i)  ||«||= 

which  completes  the  proof.  A 
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Lemma  4.4  Assume  satisfies  (H2),  (H3)  and  (H4),  and  assume  (H7)  holds, 
then  is  bounded  on  Cq. 

Proof  Let  c  be  a  constant  such  that  J^{q)  >  c^q  E  Qo  (as  guaranteed  by  (H2)). 
Assume  to  the  contrary  that  there  exists  a  point  q  E  Cq  such  that 

>  8L  (jf(5„)  -  c)  . 

Define  s  =  where  we  choose  or  small  enough  so  that  q  +  sE  Qq.  Then 

+5) - f  ( V  s)  dX  -  /'  (vjf  (5  +  As)  -  V  s)  dX 

s  tl  1^^"®  ir  0  -  f )  >  s  0  -  f )  K  -  ^)] 

This  is  positive  for  a  E  (0, 2),  thus  J^{q)  >  J^iq  +  s),  which  implies  q  +  s  eCq.  In 
addition 

^■,"(9)  -  J/'(9'  +  »)  >  J/'(9o)  -  c 

holds,  but  this  is  a  contradiction  since  q  and  9  +  s  are  in  Cq.  A 


Theorem  4.1  Assume  47^  satisfies  (H2),  (H3)  and  (H4)-  Furthermore,  assume  the 
approximate  gradient  satisfies  conditions  (H5)  and  (H6)  and  that  the  update  is  con¬ 
structed  so  that  (H7)  holds.  Then,  for  a  sufficiently  fine  discretization,  the  sensitivity 
equation  method  produces  a  sequence  of  iterates  such  that 


liminf  lli^fcll  =  0. 


k—^cx) 


(61) 


Proof  Assume  to  the  contrary  that  liminf/t^oo  |lfl'A:||  >  0  and  define  6k  such  that 


cos 


(0  \  —  (  3k  j  Sk) 

^  Iklilk.ll 


and  Wk  €.  Q  such  that 


Wk  = 


0 

1  sin(4)  (lisfcll 

sin(0fc)  =  0 


Then  {gk,  Wk)  =  0  by  construction,  and 

(si  -  V  w,)  =  -  w,)  . 
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(62) 


If  sin(0fc)  ^  0,  then  ||t<;fc||  =  1  and 


-Sfc  =  115*11  cos(^*)-jj^  +  sin(4)w*J  . 

Let  1C  denote  the  set  of  successful  iterations,  then 

„  s^jqk)  -  jr(tw)  , 

*  -  V>i(«4+l) 

for  each  k  E  )C.  Lemma  4.1  implies 


in  L*, 


Cl  Ik  II 1 


C3 


Since  is  bounded  below,  by  (H2),  the  above  condition  implies  limjk^oo.fcsx:  =  0. 
Therefore,  as  Sk  is  decreased  in  unsuccessful  iterations,  lim*_^oo  =  0.  We  now  have 
the  conditions  for  Lemma  4.2,  and 

1-  {~9ki  Sk)  _  - 

Ib^^lllki^ll " 

Thus  lim*_^oo  cos(^*)  =  1  and  lim*^oo  sin(^*)  =  0. 

Consider  the  expression 


1  —  /Ofc 


k(gfc)  -  k(gfc+i)  -  [J^{<lk)  -  J^{qk+i)) 


l-Pk< 


k(9fc)  -  kk+i)  ’ 

by  Lemma  4.3  and  the  definition  of  we  get 

1(^2  +  -  {gt  - 

{—9k,Sk)  —  I  {HkSk^Sk) 

Using  hypothesis  (H7), 

§(C2  +  i)i|s4p  -  {gt  -  VJf  fe),  s,) 

l~pk<  - - r - 

{-9k,  Sk) 

Substituting  expression  (62)  and  using  ||s*||  <  Sk,  we  get 
,  ^  2  k  +  L)Sk  -  {gk  -  {qk),9k)  -  {9k  -  ^J^iqk),  Wk)  sin(4) 

-  ||9i||cos«i 

By  Lemma  4.4  and  the  Cauchy- Schwarz  inequality,  (yiJ^{qk),‘Wk')  is  bounded  and 
we  consider  the  limit  as  k  ^  00, 

,  ,  {9k-J^(qAgk)  A9k-S^(q4 

hm  1  -  p*  =  bm  - ij— - - - — - - If. 

Ikir  Ikll 


k^co 


k^oo 
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Since  liininffc_^oo  ||fl'fc||  >  0  and  gk  is  asymptotically  consistent,  we  can  select  a  suffi¬ 
ciently  fine  discretization  such  that 

lim  1  -  pfc  <  1  -  7/2. 

K— »-00 

Hence,  pk  >  7/2  which  implies  ^*+1  >  a  contradiction.  A 


5  Duct  Design  Problem 

In  this  section,  we  use  the  duct  design  problem  to  illustrate  the  implementation  of  the 
sensitivity  equation  method.  To  begin  with,  we  will  introduce  the  discrete  approach 
for  finding  design  sensitivities  in  order  to  compare  it  with  the  sensitivity  equation 
approach. 


5.1  Discrete  Sensitivities 


To  obtain  an  algorithm  for  the  sensitivities  ^u^{q)  =  the  system  of 

nonlinear  equations  (32)  is  differentiated,  yielding 


^3+1/2  —  F i- 


3-1/2 


tpEO  _ 

^3+1/2  - 


uf,  <  Us', 

uf,  uf+1  >  Us] 
uf  <  Us  < 


(64) 


where  Fj+1/2  is  determined  by  the  scheme  used  to  compute  the  flow.  If  the  Enquist- 
Osher  scheme  was  used, 

fj+i 

h 

0 

(  d"  li-^i  Us  <  u^ , 

or  if  the  artificial  viscosity  scheme  was  used, 

=  i  (/>+' + />  -  “  -  ^“f) 

where/,-  =  /(uf,^uf), 

^  (l±a^  ‘ 


(65) 


(66) 


and 


-(  ^  A  ^ 


=  ^  (^-^A 
dq  \Adx 


u 


\Adx 


7  + 


^ 

/  dq 


u.  (67) 
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This  differentiated  scheme  can  now  be  used  to  compute 


5.2  Sensitivity  Equation 


We  now  present  the  implementation  for  the  sensitivity  equation  approach.  We  begin 
by  differentiating  the  flow  equation  (26)  with  respect  to  the  parameter  q.  Thus 


=  0 


(68) 


^m(0)  =  0  and  ^u(l)  =  0  (69) 

is  the  sensitivity  equation  for  this  problem.  Note  that  the  sensitivity  equation  is  a 
linear  equation  with  variable  coefficients  (determined  by  u).  There  has  been  little 
analysis  of  numerical  schemes  to  approximate  equations  of  this  type.  However,  for 
this  two  point  boundary  value  problem,  the  same  numerical  schemes  (Enquist-Osher 
and  artificial  viscosity)  provide  convergent  algorithms.  As  in  the  approximation  of 

to  be  the  average  sensitivity  in  the  jth  cell.  A  system  of 

nonlinear  equations  for  (^«)^(g)  =  |(^u)^(9)|  can  be  found  by  integrating 
(68)  over  each  cell. 


/  {u{xj  +  2 );  +  2))  /  (^(^i  "2)’  ~  2 )) 


3  =  1, . . . ,  iV,  where  we  assume  A  and  -^A  are  nearly  constant  over  each  cell.  As 
before,  the  terms  /  {u{xj  +  |),  §-^u{xj  +  |))  are  replaced  by  the  cell  center  values  fj 
and  /j+i.  Using  the  Enquist-Osher  scheme,  we  obtain 


and  obtain 


/i+i 

fj 

fj  +  fj+i  +  fius, 


<  Us] 
uf,uf+i  >  Us] 

uf  <  Us  <  uf^^] 

uf+l  <Us<  uf, 


(71) 


^ 

^j+1/2  -  2 


(72) 
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for  the  artificial  viscosity  scheme.  It  is  obvious  that  the  approximation  of  the  sensi¬ 
tivity  equations  depends  on  the  approximation  of  the  flow  equations.  As  described 

•  (a 

earlier,  we  use  the  notation  I  to  represent  using  scheme  N  to  approximate 

the  flow  equation  and  scheme  M  to  approximate  the  sensitivity  equation. 

5.3  Convergence  Results 

The  convergence  result  provided  in  Theorem  4.1  can  be  proved  for  the  case  when 
the  artificial  viscosity  scheme  is  used  to  approximate  the  flow  and  the  Enquist-Osher 
scheme  is  used  to  approximate  the  sensitivities  in  Algorithm  3.1.  For  this  problem, 
we  assume  (the  (Hi)  in  Theorem  4.1)  that 


^  [-^1115  ,  Qo  —  (Ajj,,  Aout)  5 

and 

=  («»)}• 

The  objective  function  given  above  is  obviously  bounded  below  (by  zero  if  all 

of  the  quadrature  weights  are  nonnegative)  satisfying  (H2).  The  hypothesis  (H3),  the 
differentiability  of 

q)  -  u{xi)y 

i=l 

on  Qo  and  hypothesis  (H4),  the  Lipschitz  continuity  of  the  derivative,  follow  from  the 
following 


Lemma  5.1  The  approximate  solution  is  differentiable  and  the  derivative  is 

Lipschitz  continuous  on  Qq. 

Proof  The  approximate  solution,  is  the  root  of  the  nonlinear  equations 

W  j)  -  +  ft  (u"-,,)  =  0,  (73) 

where  9  C°°  functions  of  their  arguments  (for  >  0).  Then  by  the 

implicit  function  theorem,  the  map 

q^u^^^{q) 


is  Lipschitz  continuously  differentiable.  A 

We  point  out  that  the  differentiability  of  the  approximate  objective  functional  is 
strongly  dependent  on  the  discretization  scheme  used  in  the  approximation.  For 
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example,  the  objective  functional  associated  with  a  Godunov  approximation  of  the 
flow  is  not  differentiable,  a  result  of  matching  a  parameter  dependent  discontinuity  on 
a  discrete  set  of  points  [4].  Finding  feasible  optimization  strategies  for  this  problem 
has  been  the  focus  of  recent  work,  see  e.g.  [4],  [19]  and  [23].  However,  for  the 
purpose  of  this  discussion,  the  artificial  viscosity  scheme  provides  a  smooth  enough 
approximate  objective  function. 

The  hypothesis  (H5)  is  guaranteed  (for  some  discretization  level)  by  the  asymptotic 
consistency  shown  below. 


Theorem  5.1  For  the  one  dimensional  Euler  equations,  the  derivative  •, 
where  the  flow  is  approximated  using  the  artificial  viscosity  approximation  and  the 
sensitivities  are  approximated  using  the  Enquist-Osher  scheme,  is  asymptotically  con¬ 
sistent  to  • 


Proof  Consider  the  norm  used  in  the  definition  of  asymptotic  consistency: 


The  first  term  on  the  right  hand  side  vanishes  since  using  the  artificial  viscosity 
scheme  for  approximating  both  the  flow  and  sensitivity  equations  leads  to  consistent 
derivatives.  The  last  two  terms  go  to  zero  as  the  approximations  Nav,  Mav  and  Meo 
are  refined,  since  the  artificial  viscosity  and  Enquist-Osher  schemes  converge  when 
used  to  approximate  the  sensitivity  equation,  jg  exact  solution  to 

the  sensitivity  equation  given  .  A 


The  hypothesis  (H6)  can  be  enforced  by  the  optimization  algorithm  by  rejecting 
steps  which  violate  this  condition  and  shrinking  the  trust-region  radius.  This  proce¬ 
dure  eventually  creates  a  step  which  satisfies  (H6),  since  the  limit  of  this  procedure 
would  produce  a  step  in  the  steepest  descent  direction. 

Finally,  (H7)  can  be  enforced  by  the  secant  update  strategy.  Therefore,  we  have 
shown  that  these  approximation  schemes  satisfy  the  conditions  of  Theorem  4.1.  Nu¬ 
merical  computations  using  these  sensitivity  schemes  are  provided  below. 
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Figure  3:  Design  Sensitivity  Approximations  Using  Enquist-Osher  Scheme 
5  A  Numerical  Results 

The  sensitivity  of  the  velocity  with  respect  to  the  Bezier  parameter,  is  presented  us- 
ing  the  numerical  schemes  described  above.  For  this  computation,  the  cross-sectional 
area  corresponds  to  an  element  of  B  (see  (35))  with  q  =  1.37125.  The  interval  [0, 1] 
is  divided  into  45  cells.  In  Figure  3,  the  sensitivity  solution  using  the  Enquist-Osher 
scheme  to  compute  both  the  flow  and  the  sensitivity  eoMbo  (compared 

with  the  closed  form  sensitivity  solution.  In  addition,  the  sensitivities  computed 
via  finite  differences  of  Enquist-Osher  solutions  using  a  finite  difference  step  size  of 
A9=  (1  X  10  q  are  also  provided.  Excellent  agreement  is  seen  for  both  of  these 
methods.  The  only  discrepancy  is  in  the  cell  to  the  left  of  the  shock,  where  numerical 
dissipation  appears  in  the  flow  solution. 

The  corresponding  design  sensitivities  which  are  computed  using  only  the  artificial 
viscosity  schemes  are  shown  in  Figure  4.  As  above,  the  agreement  is  excellent  except 
where  dissipation  errors  appear  in  the  flow  approximations.  In  this  case,  these  errors 
appear  over  more  cells  near  the  shock. 

Note  that  the  computation  of  these  sensitivities  were  performed  efficiently,  rela- 
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Figure  4:  Design  Sensitivity  Approximations  Using  Artificial  Viscosity  Scheme 

tive  to  the  cost  of  a  flow  approximation.  The  flow  approximation  requires  solving 
a  system  of  nonlinear  equations.  The  sensitivity  approximation,  on  the  other  hand, 
only  requires  solving  a  linear  system  since  the  sensitivity  appears  only  linearly  in  the 
definition  of  /  and  g.  Moreover,  if  the  Newton  method  is  used  to  solve  the  nonlin¬ 
ear  system,  then  the  linear  system  is  already  available  in  factored  form.  Therefore, 
the  sensitivities  can  be  computed  using  less  computational  time  than  required  for 
one  Newton  step.  Computational  efficiencies  such  as  this  can  be  missed  if  the  flow 
algorithm  is  simply  differentiated. 

Note  that  as  long  as  is  bounded. 


since  H  =  ul.  Thus,  one  observes  that  the  numerical  algorithms  to  compute  ei- 

■  1  d  N  f  d  \^EOi^EO 

tiler  or  are  equivalent.  This  leads  to  the  fact  that  using  the 

Enquist-Osher  scheme  to  approximate  both  the  flow  and  sensitivity  equations  pro¬ 
duces  consistent  gradients.  In  addition,  it  is  easily  seen  that  using  the  artificial 
viscosity  scheme  to  approximate  both  equations  also  produces  consistent  gradients. 
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Table  I:  A  Comparison  of  Gradients  at  the  Optimum  for  Various  Mesh  Sizes 


However,  if  the  artificial  viscosity  scheme  is  used  to  approximate  the  flow  and  the 
Enquist-Osher  scheme  is  used  to  approximate  the  sensitivity  equations,  the  gradients 
are  not  consistent  but  asymptotically  consistent. 

Numerical  results  for  this  asymptotically  consistent  case  are  provided  in  Table  I. 

6  Forebody  Simulator  Design  Problem 

We  now  describe  the  implementation  of  the  sensitivity  equation  method  for  the  fore¬ 
body  simulator  design  problem  described  in  Section  2.  As  in  the  duct  design  problem, 
we  begin  by  presenting  the  equations  which  comprise  the  discrete  sensitivity  scheme 
in  order  to  compare  and  contrast  the  two  methods.  Unlike  the  duct  problem,  we 
have  no  theoretical  convergence  results  for  the  FBS  design  problem.  However,  the 
numerical  experiments  below  show  that  the  SEM  still  converges. 

6.1  Discrete  Sensitivities 

Differentiating  the  numerical  scheme  (16)  with  respect  to  a  design  parameter,  repre¬ 
sented  by  q,  leads  to  the  following  scheme: 

[/  +  AWjA”  -  )As X 

[/  +  A«,B”  -  V,(»W  +  «'W)A,Ja,]  A^(J”  = 

[/  +  AtSr,B^  -  d-  A,  Jm]  AO” 

-  [/  -f  AtS^A^  -  -H  Ja^]  X 

Ja.  -  +  ^W)A,^ Ja,]  AQ- 


+  A<Ve(®f  -  A(Vi)A5^(J«Q*) 

+  A(V,(^®P)  -  ^¥‘>A,V,)A,(JmQ’') 

+  A«V,(»P)  -  *Wa,V,)A,^(J^(J“).  (74) 

The  equation  representing  the  boundary  conditions  are  also  differentiated.  Note 
that  the  above  sensitivity  scheme  requires  derivatives  of  the  mapping,  -^M  (denoted 
as  mesh  sensitivities)  and  the  dissipation  terms,  and  Evaluation  of 

■§^M  is  given  by  differentiating  the  scheme  which  determines  M,  see  e.g.  [20].  Other 
methods  for  approximating  have  also  been  investigated,  see  e.g.  [25].  We  see 
from  (74)  that  terms  containing  these  expressions  represent  a  significant  portion  of 
the  computational  effort,  aside  from  the  fact  that  ■§^M,  and  themselves 

need  to  be  determined. 


6.2  Sensitivity  Equation 


The  sensitivity  equation  approach  to  computing  design  sensitivities  is  presented  be¬ 
low.  To  begin  with,  we  differentiate  the  Euler  equations  and  associated  boundary 
conditions  with  respect  to  the  design  parameter  q,  which  leads  to: 


where 


dx  dy 


Fq  =  +  uQq  -t- 


Ap 

dq^ 


L  dq^"^  -^39^ 


Gq  —  -^vQ  vQq  -f 


0 
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and  where 
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■u  = 


^  (  \ 


d  pu 
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dq 


I  p  and 


d 

dq 


V  = 


\ 


d  pv 

Vj 


dq 


Ip^ 


dq 

since  /)  7^  0. 

We  are  now  free  to  apply  any  appropriate  scheme  to  solve  (75).  In  particular,  it  is 
possible  to  use  a  method  which  takes  advantage  of  the  linearity  of  the  sensitivity  equa¬ 
tion.  However,  in  this  work,  the  same  scheme  used  to  solve  the  flow  equations  is  used 
to  approximate  the  sensitivity  equations,  which  leads  to  an  ejfficient  computational 
scheme  as  in  the  discrete  approach  [2].  This  scheme  is  described  below. 

This  equation  may  now  be  transformed  to  generalized  coordinates,  so  that  the  finite 
differencing  can  be  done  more  easily.  It  makes  sense  to  use  the  same  transformation 
(which  is  equivalent  to  using  the  same  mesh)  that  was  used  in  the  solution  of  the 
Euler  equations.  Thus  the  resulting  system  is 


dK  .  dG„ 


+ 


where 


di  dq 
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where 


U  =  •  (u,  u)  and 


and 


V  =  Vq  •  (u,  u)  and 
It  can  be  shown  that 
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so  that  the  discretization  has  the  same  factored  form  as  the  Euler  equations,  thus 

[/  +  AtSfA"  -  X 

[/  +  AM,B"  -  +  4»<*>)A,J;i<]  Ai—Qr  = 

-  AtSi(F,r  -  AM,(G,)“ 

+  AiV((4'f  -  'J'f  A«V5)Aj(J';m(2)” 

+  A(V,(tf  -  4'('‘)A,V,)A,(J«(3)*  (77) 

Since  the  left  hand  side  matrices  are  the  same,  a  right  hand  side  vector  needs  to 
be  formed  for  each  design  sensitivity.  In  addition,  the  boundary  condition  type  is 
the  same  for  both  the  Euler  and  sensitivity  equations.  The  boundary  conditions  are 
determined  using  implicit  differentiation. 

Note  that  this  scheme  is  similar  to  the  discrete  sensitivity  approach.  However,  since 
the  approximation  is  applied  after  the  differentiation,  there  are  no  mesh  sensitivity 
or  dissipation  sensitivity  terms.  The  other  obvious  difference  is  that  the  boundary 
condition  on  the  parameter  dependent  boundary  is  different. 


6.3  Boundary  Conditions 

The  boundary  conditions  for  the  sensitivity  equation  (75)  are  provided  below  for  the 
case  where  the  forebody  simulator  is  described  by  a  two  parameter  Bezier  curve  (18)~ 
(20).  Extensions  to  other  forebody  descriptions  will  be  obvious.  The  appropriate 
conditions  are  obtained  by  differentiating  the  corresponding  boundary  conditions  for 
the  Euler  equations.  For  example,  at  the  inlet,  the  flow  Qin  Is  prescribed  and  will  not 
vary  as  the  forebody  parameters  q  =  (g^,  q^)  are  changed,  thus 


at  the  test  cell  inflow.  The  walls  are  treated  in  a  similar  fashion.  However,  the 
boundary  condition  at  the  forebody  simulator  surface  requires  more  attention.  This 
is  because  the  points  where  the  condition  is  evaluated  are  parameter  dependent. 

We  study  the  treatment  of  condition  (5)  in  detail.  The  normal  vector  to  the 
forebody  surface  is 


n  = 


2 

+ 


(78) 
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Thus,  the  boundary  condition  (5)  can  be  written  as 

d  d 

-u  (r„(s),  rj,(s;  q);  q)  q)  +  v  (r^(s),  rj,(s;  q)]  q)  =  0.  (79) 

The  corresponding  sensitivity  equation  boundary  condition  for  the  first  parameter, 
q^,  can  be  obtained  via  differentiation,  i.e., 

(rx(5),  Ty{s-,  q);  q)  q)  +  (rj,(s;  q),  r„(s);  q)  = 

d  d  d 

—u  (r^(s),  rj,(s;  q);  q)  9)^rj,(s;  q) 

+  u{Tx{s),Ty{s',q)]q) -^^^^Ty[s',q) 

~  9))  9) 

This  is  simply  a  nonhomogeneous  version  of  condition  (5),  namely. 


d  d 


^  ^-vJ-rT,Sr^.  (80) 


dy''dq^^^ds^^'^''dsdq^^^  dy'’ dq^^^  ds^^' 

Using  the  same  techniques,  the  boundary  conditions  corresponding  to  (6)  are; 
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The  analogous  boundary  conditions  for  q^  are  obvious. 


6.4  Numerical  Results 

The  sensitivity  equation  approach,  which  computes  design  sensitivities  for  the  two 
dimensional  Euler  equation  is  illustrated  below.  In  this  implementation,  a  right  hand 
side  vector  for  each  design  sensitivity  is  formed  along  with  the  corresponding  vector 
for  the  flow  approximations.  The  updates  for  the  flow  and  sensitivity  variables  are 
obtained  simultaneously,  exploiting  the  fact  that  the  left  hand  side  matrices  are  the 
same. 

The  design  sensitivities  with  respect  to  the  first  Bezier  parameter  q^  were  computed 
for  a  forebody  described  by  the  curve 

f  =  (x(s),y(s)),  s€[0,l]. 
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where 


xis)  =  O.OBoM  +  0ABi,3{s)  +  0.55^2, 3(s)  +  1.053,3(s), 
y{$)  =  To^o, 3(5)  +  9^51,3(5)  +  9^52,3(5)  +  FfeBa, 3(3), 

9^  =  0.1,  9^  =  0.15,  To  =  0  and  Fj  =  0.2.  This  curve  is  twice  as  long  in  the  s-direction 
as  the  admissible  forebody  simulators  given  in  B  (see  (18)).  Under  a  uniform  inlet 
flow  profile  described  by  the  inlet  Mach  number,  Ma  =  2.0,  the  approximate  flow 
variables  and  sensitivities  are  computed  on  a  43  x  49  mesh.  The  sensitivity  of  the 
x-component  of  momentum  with  respect  to  the  Bezier  parameter  9^,  computed  using 
the  sensitivity  equation  approach  and  the  finite  difference  approach  (for  4  different 
step  sizes)  are  plotted  along  the  outflow  plane  in  Figure  5.  The  corresponding  plots 
for  the  Energy  sensitivity  are  provided  in  Figure  6.  Observe  that  the  step  size  of 
0.00001  produces  noisy  sensitivity  values  close  to  the  forebody  (presumably  due  to 
round-off  errors).  A  larger  step  size  of  0.01  gives  the  best  results  (when  compared 
to  the  sensitivity  equation  approach)  near  the  shock  location.  The  best  qualitative 
behavior  appears  when  the  step  size  is  0.001.  These  figures  demonstrate  the  difficulty 
of  obtaining  a  satisfactory  step  size  at  all  resolution  levels  in  the  flow  domain. 

A  model  forebody  simulator  design  problem  is  discussed  below.  To  begin  with, 
we  seek  the  optimum  value  of  the  inlet  Mach  number  and  two  Bezier  parameters 
(  (9S9^)»  describing  a  shortened  forebody  simulator  in  the  admissible  set  B)  which 
minimize  the  approximate  cost  functional  (given  in  equation  (25)).  The  flow 
data  Q  to  be  matched  is  given  by  the  flow  corresponding  to  the  forebody  shape 
F  described  above.  We  point  out  that  the  artificial  dissipation  in  the  flow  solver 
produces  a  “smearing”  effect  on  the  flow  variables.  Therefore,  based  on  the  results 
for  the  duct  design  problem,  we  expect  a  sufficiently  smooth  approximate  cost  func¬ 
tional.  Furthermore,  the  comparison  of  the  sensitivities  in  Figures  5  and  6  lead  us  to 
believe  that  the  sensitivity  equation  approach  may  produce  asymptotically  consistent 
derivatives. 

The  sensitivity  equation  method  was  applied  to  the  FBS  design  problem  with 
initial  values  of  the  parameters;  Ma  =  2.0,  9^  =  0.10  and  9^  =  0.15.  These  parameters 
correspond  to  those  used  to  generate  Q  (even  though  that  forebody  is  longer).  We 
present  the  iteration  history  in  Table  II.  Observe  that  there  is  a  drastic  reduction 
in  the  approximate  cost  functional  in  the  first  three  iterations.  The  iteration  history 
for  the  x-component  of  momentum  is  given  in  Figure  7.  Note  that  the  front  end  of 
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Figure  5:  Comparison  of  Sensitivities  at  Outflow:  x-Component  of  Momentum 
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Figure  6:  Comparison  of  Sensitivities  at  Outflow:  Energy 
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Table  II:  Shortened  Forebody  Optimization 


Iteration 

Ma 

Q2 

Cost  Functional 

Gradient 

0 

2.00000 

0.10000 

0.15000 

3.2339 

27.1283 

1 

2.00108 

0.14608 

0.17177 

1.6000 

11.6285 

2 

2.01054 

0.26846 

0.14152 

0.3332 

3.7955 

3 

2.00897 

0.30765 

0.13671 

0.2334 

0.4621 

4 

2.01027 

0.30139 

0.14007 

0.2306 

0.5963 

5 

2.01307 

0.29367 

0.14737 

0.2289 

0.6861 

6 

2.01670 

0.28891 

0.15564 

0.2271 

0.5009 

7 

2.01900 

0.29011 

0.15921 

0.2249 

0.1513 

8 

2.01940 

0.29278 

0.15821 

0.2237 

0.0576 

9 

2.01936 

0.29420 

0.15669 

0.2232 

0.0571 

10 

2.01952 

0.29439 

0.15604 

0.2230 

0.0275 

11 

2.01994 

0.29417 

0.15603 

0.2229 

0.0173 

12 

2.02006 

0.29415 

0.15609 

0.2229 

0.0153 

the  forebody  simulator  becomes  more  blunt  during  the  first  two  iterations  in  which 
a  stagnation  region  is  set  up  in  front  of  the  FBS.  This  has  the  effect  of  moving  the 
shock  forward,  which  comes  close  to  the  shock  location  created  by  the  long  forebody. 
The  remaining  iterations  are  used  to  “fine  tune”  the  solution  near  the  FBS.  The 
comparison  of  the  optimal  forebody  simulator  to  the  flow  generated  by  the  long 
forebody  is  displayed  in  Figure  8.  Notice  that  the  shock  location  is  the  same  in  both 
flows. 

In  the  optimization  above,  the  initial  Hessian  was  computed  using  forward  dif¬ 
ferences.  This  adds  some  initial  expense  in  the  hope  for  fewer  iterations.  However, 
without  this  technique,  using  the  identity  matrix  as  the  initial  Hessian,  the  iteration 
converged  in  fifteen  iterations.  Therefore,  neither  technique  showed  an  advantage. 

6.5  Conclusions 

While  no  rigorous  proof  of  asymptotically  consistent  gradients  has  been  shown  for 
Euler  equations,  numerical  evidence  in  [3]  suggests  that  the  gradients  may  indeed 
be  asymptotically  consistent.  Similar  numerical  evidence  exists  for  finite  element 
approximations  of  the  Navier-Stokes  equations  [5]. 
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Iteration  4  Iteration  6  Iteration  8  Iteration  12 

Figure  7:  Sensitivity  Equation  Method  Iterations 
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Final  Converged  Solution  Original  Long  Forebody  (data  to  be  matched) 

Figure  8:  Comparison  of  Optimal  Solution  with  Data 


References 

[1]  R.  Beam  and  R.  F.  Warming,  Journal  of  Computational  Physics  22,  87 
(1976). 

[2]  J.  Borggaard,  J.  Burns,  E.  Cliff,  and  M.  Gunzburger,  Sensitivity  Cal¬ 
culations  for  a  2D,  Inviscid,  Supersonic  Forebody  Problem,  in  Identification  and 
Control  of  Systems  Governed  by  Partial  Differential  Equations,  edited  by  H.  T. 
Banks,  R.  Fabiano,  and  K.  Ito,  pp.  14-24,  Philadelphia,  PA,  1993,  SIAM 
Publications. 

[3]  J.  T.  Borggaard,  The  Sensitivity  Equation  Method  for  Optimal  Design,  PhD 
thesis,  Virginia  Tech,  Blacksburg,  VA,  1994. 

[4]  J.  T.  Borggaard,  On  the  Presence  of  Shocks  in  Domain  Optimization  of  Euler 
Flows,  in  Flow  Control,  edited  by  M.  Gunzburger,  volume  68  of  Proceedings 
of  the  IMA,  Springer- Verlag,  1995. 

[5]  J.  V.  Burkardt,  Sensitivity  Analyses  and  Computational  Shape  Optimization 
for  Incompressible  Flows,  PhD  thesis,  Virginia  Tech,  Blacksburg,  VA,  1995. 

[6]  J.  A.  Burns,  B.  B.  King,  and  Y.-R.  Oh,  A  Computational  Approach  to  Sen¬ 
sor/Actuator  Location  for  Feedback  Control  of  Fluid  Flow  Systems,  in  Sensing,, 
Actuation,  and  Control  in  Aeropropulsion,  edited  by  J.  D.  PaduanO,  volume 
2494  of  Proceedings  of  the  International  Society  for  Optical  Engineering,  pp.  60— 
69,  1995. 

[7]  R.  G.  Carter,  Numerical  Optimization  in  Hilbert  Space  Using  Inexact  Func¬ 
tion  and  Gradient  Evaluations,  Technical  Report  89-45,  ICASE,  1989. 

[8]  R.  G.  Carter,  SIAM  Journal  of  Numerical  Analysis  2S,2I>\  (1991). 

[9]  G.  K.  Cooper  and  J.  R.  SirBAUGH,  PARC  Code:  Theory  and  Usage,  Tech¬ 
nical  Report  AEDC-TR-89-15,  Arnold  Engineering  Development  Center,  Arnold 
AFB,  TN,  1989. 

[10]  J.  E.  Dennis  Jr.,  D.  M.  Gay,  and  R.  E.  Welsch,  Transactions  on  Mathe¬ 
matical  Software  7,  348  (1981). 

[11]  J.  E.  Dennis  Jr.  and  R.  B.  Schnabel,  Numerical  Methods  for  Unconstrained 
Optimization  and  Nonlinear  Equations,  Prentice  Hall,  Englewood  Cliffs,  New 
Jersey,  1983. 

[12]  G.  FaRIN,  Curves  and  Surfaces  for  Computer  Aided  Geometric  Design:  A  Prac¬ 
tical  Guide,  Academic  Press,  San  Diego,  CA,  1988. 

[13]  P.  D.  Frank  and  G.  R.  Shubin,  Journal  of  Computational  Physics  98,  74 
(1992). 


37 


[14]  E.  J.  Haug,  K.  K.  Choi,  and  V.  Komkov,  Design  Sensitivity  Analysis  of 
Structural  Systems,  volume  177  of  Mathematics  in  Science  and  Engineering,  Aca¬ 
demic  Press,  Orlando,  FL,  1986. 

[15]  D.  Huddleston,  Development  of  a  Free- Jet  Forebody  Simulator  Design  Opti¬ 
mization  Method,  Technical  Report  AEDC-TR-90-22,  Arnold  Engineering  De¬ 
velopment  Center,  Arnold  AFB,  TN,  1990. 

[16]  A.  Jameson,  Journal  of  Scientific  Computing  3,  233  (1988). 

[17]  A.  Jameson  and  J.  Reuther,  Control  theory  based  airfoil  design  using  the 
Euler  equations,  in  Proceedings  from  the  5th  AIAA/USAF/NASA/ISSMO  Sym¬ 
posium  on  Multidisciplinary  Analysis  and  Optimization,  pp.  206-222,  1994. 

[18]  C.-Y.  JOH,  B.  Grossman,  and  R.  T.  Haftka,  Engineering  Optimization  2\, 
1  (1993). 

[19]  R.  NarDUCCI,  B.  Grossman,  and  R.  Haftka,  Inverse  Problems  in  Engineer¬ 
ing  2,  49  (1995). 

[20]  N.  Pagaldipti  and  A.  Chattopadhyay,  A  Discrete  Semi- Analytical  Proce¬ 
dure  for  Aerodynamic  Sensitivity  Analysis  Including  Grid  Sensitivity,  in  Proceed¬ 
ings  from  the  5th  AIAA/USAF/NASA/ISSMO  Symposium  on  Multidisciplinary 
Analysis  and  Optimization,  pp.  161-169,  1994. 

[21]  O.  Pironneau,  Journal  of  Fluid  Mechanics  59,  117  (1973). 

[22]  O.  Pironneau,  Journal  of  Fluid  Mechanics  60,  97  (1974). 

[23]  A.  Shenoy  and  E.  M.  CLIFF,  An  optimal  control  formulation  for  a  flow  match¬ 
ing  problem,  in  Proceedings  from  the  5th  AIAA/USAF/NASA/ISSMO  Sympo¬ 
sium  on  Multidisciplinary  Analysis  and  Optimization,  pp.  520-528,  1994. 

[24]  A.  C.  Taylor  III,  G.  W.  Hou,  and  V.  M.  Korivi,  A  Methodology  for  Deter¬ 
mining  Aerodynamic  Sensitivity  Derivatives  with  Respect  to  Variation  of  Geo¬ 
metric  Shape,  in  Proceedings  of  the  AIAA/ASME/ASCE/AHS/ASC  32nd  Struc¬ 
tures,  Structural  Dynamics,  and  Materials  Conference,  Baltimore,  MD,  1991, 
AIAA  paper  91-1101. 

[25]  A.  C.  Taylor  III,  G.  W.  Hou,  and  V.  M.  Korivi,  Sensitivity  Analysis,  Ap¬ 
proximate  Analysis  and  Design  Optimization  for  Internal  and  External  Viscous 
Flows,  in  Proceedings  of  AIAA  Aircraft  Design  Systems  and  Operations  Meeting, 
1991,  AIAA  paper  91-3083. 

[26]  J.  Thompson,  Z.  Warsi,  and  C.  Mastin,  Numerical  Grid  Generation:  foun¬ 
dations  and  applications,  North-Holland  Publishing  Company,  New  York,  1985. 


38 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No,  0704-0188 


Pu  Wic  wporUnj  burden  for  thif  collection  information  it  estimated  to  average  1  hour  per  response,  includingthc  time  for  reviewing  instructions,  searching  existing  data  sources, 
gaw^ngand  maintain  mg  the  data  needed,  and  compicbngand  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this 

^  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports.  1215  Jefferson 
Davis  Highway,  Suite  1204,  Arlington.  VA  22202-4302,  and  to  the  Office  of  Management  and  Budget.  Papenwrk  Reduction  Project  (0704-0188),  Washington,  DC  20503 


1.  AGENCY  USE  OHlY(Leave  blank)i2.  REPORT  DATE 

I  Jnne  1996 


3.  REPORT  TYPE  AND  DATES  COVERED 
Contractor  Report 


4.  TITLE  AND  SUBTITLE 

A  PDE  SENSITIVITY  EQUATION  METHOD 

FOR  OPTIMAL  AERODYNAMIC  DESIGN 

5.  FUNDING  NUMBERS 

C  NASl-19480 

WU  505-90-52-01 

6.  AUTHOR(S) 

Jeff  Borggaaxd  and  John  Burns 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  AOORESS(ES) 

Institute  for  Computer  Applications  in  Science  and  Engineering 

Mail  Stop  132C,  NASA  Langley  Research  Center 

Hampton,  VA  23681-0001 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

ICASE  Report  No.  96-44 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

National  Aeronautics  and  Space  Administration 

Langley  Research  Center 

Hampton,  VA  23681-0001 

10.  SPONSORING/MONITORING 

AGENCY  REPORT  NUMBER 

NASA  CR-198349 

ICASE  Report  No.  96-44 

11.  SUPPLEMENTARY  NOTES 

Langley  Technical  Monitor:  Dennis  M.  Bushnell 

Final  ^port 

Submitted  to  the  Journal  of  Computational  Physics. 

12*.  OISTRIBUTION/AVAILABILITY  STATEMENT 

Undassified-Unlimited 

12b.  DISTRIBUTION  CODE 

Subject  Category  64 

13.  ABSTRACT  (Mejdmum  200  words)  | 

The  use  of  gradient  based  optimization  algorithms  in  inverse  design  is  well  established  as  a  practical  approach 
to  aerodynamic  design.  A  typic^  procedure  uses  a  simulation  scheme  to  evaluate  the  objective  function  (from  the 
approximate  states)  and  its  gradient,  then  passes  this  information  to  an  optimization  algorithm.  Once  the  simulation 
scheme  (CFD  flow  solver)  has  been  selected  and  used  to  provide  approximate  function  evaluations,  there  are  several 
possible  approaches  to  the  problem  of  computing  gradients.  One  popular  method  is  to  differentiate  the  simulation 
scheme  and  compute  design  sensitmties  that  are  then  used  to  obtain  gradients.  Although  this  black-box  approach 
has  many  advantages  in  shape  optimization  problems,  one  must  compute  mesh  sensitivities  in  order  to  compute 
the  design  sensitivity.  In  this  paper,  we  present  an  alternative  approach  using  the  PDE  sensitivity  equation  to 
develop  algorithms  for  computing  gr^ents.  This  approach  has  the  advantage  that  mesh  sensitivities  need  not  be 
computed.  Moreover,  when  it  is  possible  to  use  the  CFD  scheme  for  both  the  forward  problem  and  the  sensitivity 
equation,  then  there  are  computational  advantages.  An  apparent  disadvantage  of  this  approach  is  that  it  does  not 
always  produce  consistent  derivatives.  However,  for  a  proper  combination  of  discretization  schemes,  one  can  show 
asymptotic  consistency  under  mesh  reflnement,  which  is  often  sufficient  to  guarantee  convergence  of  the  optimal 
dengn  algorithm.  In  particular,  we  show  that  when  asymptotically  consistent  schemes  are  combined  with  a  trust- 
region  optimization  algorithm,  the  resulting  optimal  design  method  converges.  We  denote  this  approach  as  the 
sensitivity^  equation  method.  The  sensitivity  equation  method  is  presented,  convergence  results  are  given  and  the 
approach  is  illustrated  on  two  optimal  design  problems  involving  shocks. 


14.  SUBJECT  TERMS 

Optimal  Design;  Sensitivity  Equations;  Trust  Region  Methods 


17.  SECURITY  CLASSIFICATION 
OF  REPORT 
Unclassified 


7540-0l-2a0-SSd0 


15.  NUMBER  OF  PAGES 
40 


16.  PRICE  CODE 

_  A03 


16.  SECURITY  CLASSIFICATION  19.  SECURITY  CLASSIFICATION  20.  LIMITATION 
OF  THIS  PAGE  OF  ABSTRACT  OF  ABSTRACT 

Unclassified 


^^^^^^^^^^^sSndarc^oriT^3!7Rev^^ 

Ptescribed  by  ANSI  Std.  Z39-18 
298-102 


